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I  13.  ABSTRACT  200 


A  successful  nonlinear  partial  difr_rential  equation  based  approach  to  restoration 
was  carried  out,  ENO  least  squares,  shock  filters,  feature  detectors  and  total 
variation  based  deconvolut i  'u  techniques  were  combined.  Also  rigorous  morphological 
methods  and  wavelet  analysi.;  were  developed  and  used  to  restore  noisy,  blurry 
images . 
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I.  Introduction. 


A  nonlinear  partial  differential  beised  approach  to  some  of  the  basic  problems  of 
image  processing  was  initiated.  Problems  of  noise  removal,  enhzincement,  and  ap¬ 
proximation  of  restored  noisy,  blurry  images  were  attacked  using  this  new  approach. 
We  have  overcome  many  of  the  difficulties  experienced  by  standard  techniques  such 
as  spurious  oscillations  (ringing)  and/or  excessive  blurring  of  edges.  Here,  we  used 
essentially  nonoscillatory  (ENO)  least  squares  approximation,  (see  e.g.  [3]),  to¬ 
gether  with  feature  detectors,  shock  filters  [4]  and  total  variation  [5,  6,  21]  based 
deconvolution  to  produce  a  state-of-the-zu’t  enhancement  technique.  Moreover  the 
ENO  preprocessing  was  used  together  with  a  flame  filter  based  on  Hamilton- Jacobi 
type  ideas.  This  was  generalized  to  the  unique  class  of  morphological  requirements 

[19]  satisfying  operators  using  numerical  techniques  simileir  to  those  originating  in 
[2].  Additionailly,  an  amisotropic  diffusion  process  originated  by  Malik  and  Perona 

[20]  was  turned  into  an  anisotropic  shock/rarefaction  filter.  The  resulting  evolution 
equation  both  enhances  and  denoises. 

A  major  breakthrough  came  in  our  total  variation  based  restoration  of  noisy, 
blurred  images.  The  total  variation  of  the  image  is  minimized  subject  to  constraints 
involving  the  point  spread  function  of  the  blurring  process  and  the  statistics  of  the 
noise.  The  solution  is  obtained  using  Euler- Lagrange  equation  with  artificieil  time 
evolution,  and  Lagrange  multipliers  enforce  the  constraints.  This  amounts  to  solving 
an  interesting  time  dependent  partial  differential  equation  on  a  manifold  determined 
by  the  constraints.  As  t  increases  the  image  is  restored.  The  numerical  algorithm 
is  simple  to  implement  and  appears  to  be  nonoscillatory  (minimal  ringing)  and 
noninvasive  (recovers  sharp  edges). 

II.  Restoration  Algorithms. 

Our  first  restoration  algorithm  involved  additive  noise; 

We  solve  the  following  problem.  Let 


(2.1) 


«o(a:,  y)  =  (Au)(i,  y)  -f  n(i,  y) 


where  A  is  a  linear  integral  operator  and  n  is  additive  noise.  Also  uo{x,y)  is  the 
observed  intensity  function  while  u{x,y)  is  the  image  to  be  restored.  Both  are 
defined  on  a  region  f2  in  /2^.  For 

NTTS 


A  may  be  a  convolution  type  integral  operator  in  which  case  we  write: 


(2.2)  (Au)(x,  y)  =  (A:  *  u)(x,  y). 

Examples  we  shall  experiment  with  include 


lAB  n 

i  Un«i‘iao3Ji#84  Q 
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Motion  blur: 


(2.3a) 


k{x,y)  =  -  if  — 

^  a  2  ~  ~  2 

y)  =  0  elsewhere. 


Diffraction  limited  blur: 

(2.3b)  k{x,y)  =  a  sine  ax 

where  sine  x  =  and  a  is  related  to  the  aperture  size  and  bandwidth  limitation 
of  the  transmission  system. 

Defocus  blur: 


(2.3c) 


Ji(aWwf+^) 

k{u:^,ujy)  =  - y  - 

aJu^l+wl 


where  *  denotes  the  Fourier  transform  Ji(x)  is  the  Bessel  function:  and  a  >  0 
measures  the  degree  of  blurring  (a  =  0  is  just  the  identity  map). 

Gaussian  blur: 


(2.3d) 


k{x,  y)  =  (Ana)  ^  exp 


(x^  +  y^) 

4a 


which  arises  e.g.,  in  atmospheric  turbulence. 

All  these  kernels  result  in  severely  blurred  images  whose  restoration  is  an  ill-posed 
problem.  Inverting  the  equation 

(2.4)  k  *u  =  Uo 

leads  to  the  formula 

u  =  J^-\k)-'^J^Uo 

where  !F  is  the  Fourier  transform,  .F~‘  its  inverse  and  k  is  the  transform  of  the 
kernel.  The  ill-posedness  comes  from  the  zeros  of  k(^,  rf)  for  rj  in  bounded  regions 
of  and  for  —*  oo.  Even  in  the  absence  of  noise,  quantization  error  leads 

to  severe  ringing  in  this  process. 
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The  ill-posedness  of  this  procedure  has  long  been  well  known.  Conventioneil 
variational  approaches  to  this  problem  involve  a  least  squares  fit  because  this 
leads  to  lineeir  equations.  The  first  attempt  along  these  lines  was  made  by  Phillips 
[11]  and  later  refined  by  Twomey  [12,13]  in  the  one  dimensional  case.  In  our  two 
dimensional  continuous  framework  their  constrained  minimization  problem  becomes 


(2.5a) 


Minimize 


+  Uyyfdxdy 


subject  to  constraints  involving  the  meein 

(2.5b)  f  udxdy  =  I  uodxdy 

j\i  Ju 

and  standard  deviation 

(2.5c)  /  {Au  —  Ufii^dxdy  —  . 

Ju 


The  resulting  linear  system  is  easy  to  solve  using  modern  numerical  linear  algebra. 
However,  the  results  are  rather  disappointing,  see  section  VI  below. 

We  instead  minimize  the  variation  of  the  image,  which  is  a  direct  measure  of  how 
oscillatory  it  is.  The  space  of  functions  of  bounded  variation  (BV)  is  appropriate 
for  discontinuous  functions.  This  is  well  known  in  the  field  of  shock  calculations, 
e.g.,  [3],  and  references  therein. 

Thus,  our  constrained  minimization  problem  is 


(2.6a) 


Minimize 


L-f 


I  +  u'^dxdy 


subject  to  constraints  involving  uq. 

In  our  work  so  far  we  have  taken  the  same  two  constraints  (2.5b)  8md  (2.5c). 

Again  (2.5b)  indicates  that  the  white  noise  n(x,y)  in  (2.1)  is  of  zero  mean,  and 
(2.5c)  uses  a  priori  information  that  the  standard  deviation  of  the  noise  n(x,  y)  is 

<7. 

In  this  case  we  arrive  at  the  Euler- Lagrange  equations 


(2.7a) 


—  ^(A*Au  —  A*u)  in  Q 
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with 


(2.7b) 


=  0  on  Q. 
an 


Here  .4*  is  just  the  adjoint  integral  operator. 

The  constant  A  is  a  Lagrange  multiplier  chosen  so  that  the  constraint  (2.5b)  is 
satisfied.  The  first  constraint  (2.5a)  will  be  shown  in  Remark  4  below  to  be  satisfied 
if; 


(2.8) 


y^.4u = y*  “  J  ^  J  ^ 


for  each  u.  This  is  true  up  to  norm?dization  for  A  a  convolution.  Thus  we  assume 
(2.8)  for  simplicity  only. 

We  shall  use  the  gradient-projection  method  of  Rosen  [15]  which,  in  this  case, 
becomes  the  interesting  “constrained”  partial  differential  equation 


(2.9a) 


Ut  = 


—  XA*{Au  —  uq) 


for  <  >  0,  (x,y)  e  fl 


(2.9b) 


P  =  0  on  afi 
On 


aind  u(x,t/,0)  given  so  that  (2.5b,c)  eire  satisfied.  If  (2.5b)  is  satisfied  initially,  e.g. 
u(i,y,0)  =  Uo(x,y),  then,  by  the  conservation  form  of  (2.9a,b)  eind  by  (2.8),  it  is 
always  satisfied.  Satisfying  (2.5c)  can  be  done  through  a  process  described  in  the 
next  section. 

The  projection  step  in  the  gradient-projection  method  just  amounts  to  updating 
A(t)  so  that  (2.5c)  remains  true  in  time.  One  merely  differentiates  (2.5c)  with 
respect  to  time,  replaces  Ut  by  the  right  side  of  (2.9)  and  then  chooses  A(t)  so 
that  this  is  set  to  be  zero.  This  amounts  to  multiplying  the  right  side  of  (2.9)  by 
A*(Au  —  uo),  integrating  over  f2,  setting  the  result  equal  to  zero,  and  solving  for 
A(t). 
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We  have 


»  ( (  v/4+^)  )  )  ‘ 

J^(A*(Au  -  uo)ydxdy 

We  thus  have  a  dynamical  procedure  for  restoring  the  image.  As  t  — »  oo  the 
steady  state  solution  is  the  desired  restoration. 

Remark  (1.1).  Equation  (2.9a)  can  be  written  as 
(2.11)  ut  =  IC{u)  -  \A*iAu  -  uo) 

where  AC(u)  is  the  curvature  of  the  level  set  u  =  constant  at  each  point.  This  part 
of  the  operator  can  be  viewed  as  moving  each  level  set  normal  to  itself  with  velocity 
equal  to  its  curvature,  divided  by  the  magnitude  of  the  gradient.  The  constraint 
term  just  acts  to  project  the  motion  back  so  that  (2.5c)  is  satisfied. 

Remark  (1.2).  The  method  is  quite  general  as  regards  nature  and  number  of 
constraints.  Generally  one  has  many  Lagrange  multipliers  and  one  must  invert  a 
Gram  type  of  matrix  to  update  their  values.  In  particular,  one  may  localize  the 
constraints  over  space. 

Remark  (1.3).  The  method  is  also  general  as  regards  the  nature  of  the  noise.  We 
have  recently  experimented  successfully  with  both  multiplicative  and  speckle  noise. 
Again,  we  demonstrate  the  former  in  section  VI  and  discuss  the  algorithm  in  section 
IV. 

III.  Numerical  Approximations. 

The  easiest  way  to  construct  an  efficient  numerical  approximation  to  (2.7)  is  to 
approximate  the  variational  problem  (2.6)  with  constraints  (2.5b,c)  directly. 

We  do  this  as  follows:  Let  Xi  =  iAi,  yj  =  j  Ax,  i,j  =  0, 1, . . .  ,  iV,  with  NAx  =  1, 
and  also  Uij  ~  u{iAx,jAx).  Our  discrete  variational  problem  is 


N-i  _ 

(3.1a)  Minimize 

«,j=0 

subject  to  the  constraints: 


(3.1b) 


N  N 

UijAx  =  (uo)i>Aa: 

i,j=0  i,j=0 
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(3.1c) 


N 


Mu  -  uof  =  Yi  l('4u)iy  -  (uo)oP(Ax)  = 

»,j=0 


Here  A^Ujj  =  ±(ui±ij  —  similarly  for  ?ind  e  >  0,  of  the  order  of 

round-off  error,  is  used  to  avoid  division  by  zero. 

The  discrete  Euler- Lagrange  equations  come  just  by  differentiating  (3.1)  with 
respect  to  Uij  and  using  a  Lagrange  multiplier.  We  arrive  at: 


(3.2a) 


0  =  A1 


■^(A5.u.j)2-|-(A![.u.;)2  -f-c 

_ ^>0 _ 

Y(A- tx.,)2 -KA!' u,,)2  + 
-\[{A*Au)ij-{A*u)ij] 


+  Al 


with  boundary  conditions 


(3.2b)  A^uoj  =  0  =  A^Uiv-i.j,  j=0, 1,...,A^ 

(3.2c)  A!(.u,o  =  0  =  A^u,,/v_i,  i  =  0,l,...,N 


We  approximate  this  by  the  gradient-projection  method: 


(3.3) 


+  c 


A1 


+  A?L 


y/{A%Uijy  +  (A!|.Uij)2  -I- 


A^ 


y^(A5.«,j)2  -h(A!}.Ui>)2 
\[(A*Au)^j-iA-uoh]. 


Typical  values  c  =  ^  range  from  .3  to  3. 

Remark  (1.4).  We  also  used  a  more  isotropic  approximation  of  the  discrete  vari- 
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ation;  i.e.  minimize 


(3.1c) 

*J=1 

r»  — 1 

t  i 

where 

(3.1e) 

~  V2 

(3.1f) 

A,'),,  _  ~  «i> 

yielding  simple  expressions  corresponding  to  (3.2a)  and  (3.3).  The  results  are 
slightly  more  pleasing. 

Rather  than  use  a  discrete  version  of  (2.10)  we  impose  (3.1c)  directly  on  in 
(3.3).  This  amounts  to  solving  a  quadratic  equation  for  A"  as  follows 

(3.4a)  -  uof  = 


This  means,  among  other  things,  that  the  constraint  (3.1c)  need  not  be  satisfied 
by  the  initial  guess  In  our  experiments  we  usually  took  =  (uo),j.  Also 

the  resulting  quadratic  equation  for  A  derived  below  has  (in  general)  two  complex 
roots: 


(3.5) .  A±  =  Cl  ±  (sign  ci)v^ 

If  C2  <  0,  we  just  set  A  =  cj.  If  C2  >  0  we  take  A  =  A_.  We  have  found  in 
our  experiments  that  this  procedure  eventually  drives  us  to  the  manifold  (3.1c)  for 
<T  >  0.  Constraint  (3.1b)  is  automatically  satisfied  by  the  procedure  (3.3)  when  the 
initial  data  satisfies  the  constraint. 

To  simplify  the  solution  procedure  we  rewrite  (3.3)  as 

(3.6)  -  A"g.';. 

and  define  {p"  }  =  p",  {9"  )  =  9”. 


7 


Then  (3.4a)  becomes 


(3.7) 


-  uolp  =  =  Pp"  -  AMg"  -  uof 

=  A2pg"||2  -  2A(>lg«,^p”  -  «o)  +  IMp"  - 


(dropping  the  superscript  n). 


A±  =  (ll>igll)  ^  sign(>lg,  Ap  -  uq) 


(3.8) 


l(Ag,Ap-uo)| 


±  sign(Ag,  Ap  -  uo)\/(Ag,Ap  -  uq)^  -  ||Ag||2(||Ap  -  uqIP  - 


The  definition  of  A_  above  approaches  that  obtained  by  the  gradient-projection 
method  in  (2.10).  Precisely  we  may  write 


A_  = 


(3.9) 


1  (||Ap-uo|p  -g^) 

2  (Aq,Ap~uo) 

,  ^nW(Pp-tiolP-cr-)^ 

V  l(-4g,Ap- uo)^! 


Using  (3.6),  (3.3),  and  Taylor’s  theorem  we  have  (in  smooth  regions) 
(3.10) 

r  r 

d  I  Uj 
dr 


(Ap  -  uo)tj  ~  (Au),j  -  {uo)ij  -I-  At 


d 

+ 


'At  +  ul 


+  0(Ax) 


So  if  (3.10)  is  satisfied  at  time  level  n  we  have 


(3.11) 


Also 


||Ap  —  uoll^  —  cr^  =  2At(Au"  —  uo,A[ 


dx 


Ux  O  Uy 

yjul  +  Uy  \/^z  + 


]) 


+  0((At)2-|-(Ax)2) 


(3.12) 


(Aq)ij  =  (AA*  ■  (Au  -  uo))ij. 
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So 


(3.13) 


{Aq,Ap  -  uo)  =  {A*{Au^  -  uo),A*{Au''  -  uo)) 
+  0{At). 


Finally,  from  (2.10),  (3.11),  (3.13)  and  (3.9)  we  have 

A" 

(3.14)  ^  =  A(nA0  +  O(A<  +  Ax) 


which  is  the  desired  result. 

IV.  Multiplicative  Noise. 

This  time  we  shall  solve  the  following  problem.  Let 

(4.1)  UQ(x,y)  =  nAu. 

where  A  is  as  in  (2.1)  here  n  has  mean  one  and  standard  deviation  . 

We  shall  use  the  same  cleiss  of  convolution  kernels  as  in  the  additive  noise  case. 
We  shall  again  minimize  the  variation  of  the  image,  as  in  (2.6c),  this  time  subject 
to  constraints  of  the  following  type; 


(4.2a) 


and 


(involving  the  mean) 


(4.2b) 

1  f  {Au  ~  uq)^ 

2/0  rUo)2 


dxdy  = 


(involving  the  standard  deviation) 


The  resulting  gradient-projection  algorithm  involves  two  Lagrange  multipliers. 
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We  follow  the  gradient-projection  technique,  solving  for  both  parameters  with 
very  successful  results,  as  we  shall  see  in  section  VI  below. 


5.  Theoretical  just’acation. 

This  work  is  joint  with  Professor  Pierre- Louis  Lions  We  begin  with  einalysis  of 
the  minimization  problem  (2.6a)  subject  to  linear  constreiints. 

Vv’e  first  need  to  formulate  the  problem  in  BV{Q.)  -  we  assume  that  f2  is  smooth 
enough  so  that  extensions  to  BV'(R^)  are  possible,  a  fact  that  is  true  as  long  as  Q 
is  smooth  or  if  fl  is  a  rectangle. . . .  We  denote  by  lL)u|  the  seminorm  on  BV 

that  coincides  with  yjul  +  u^dxdy  when  u  is  smooth  (or  . . . ).  Then,  we 
consider  the  following  minimization  problem: 

(5.1)  Inf  {  f  \Du\/u  e  J3V'(Q),  /  |j4u  —  uol^  =  cr^,  /  (u  -  uo)  =  0). 

Ju  Jii  Jn 


Proposition  5.1:  The  minimization  problem  (5.1)  admits  at  least  one  solution  if 
A  is  compact  from  L^{U)  into  L^{Q). 

Proof  We  first  aissume  that  because  of  (2.8)  A1  =  1  and  thus  by  standard  functional 
analysis  arguments,  |Dul  -i-  ||.4u||£,2(n)  is  a  norm  on  BV(Q)  which  is  equivalent 
to  the  usual  norm. 

Then,  this  implies  that  minimizing  sequences  of  (5.1)  are  bounded  in  B\''{^) 
and  thus  in  L^{U)  by  Sobolev  imbedding.  Therefore,  if  {u„}  denotes  am  arbitrary 
minimizing  sequence  of  (5.1),  we  may  asume  without  loss  of  generality  that  {un) 
converges  weakly  to  some  u  in  BV{n)  (weak-*  convergence)  and  in  L^(n).  There¬ 
fore,  we  recover  at  the  limit  (2.5b)  and  (2  5c)  follows  from  the  assumption  that  A 
is  compact  so  that  .4u„  converges  to  Au  in  L^{Vt).  We  conclude  then  easily  that  u 
is  a  minimum  since  wc  have  by  weak  convergence 


/  |Du|<  —  /  |Du„|. 
Jn  «  Jn 


□ 


Remark  5.1:  Notice  that  the  assumption  made  in  Proposition  (5.1)  on  A  excludes 
the  obviously  interesting  case  when  A  is  the  identity  operator,  in  which  case  (5.1) 
turns  out  to  be  a  highly  non-trivial  minimization  problem  related  to  isoperimetric 
inequalities  and  geometrical  problems.  Any  discussion  of  this  would  be  too  technical 
for  the  main  purpose  of  this  report. 


We  next  study  equations  of  the  form  (2.9a)  namely 


^  0 

(5.2a)  ut  =  p — {a,(Vu)}  —  XA*{Au  —  uo)  for  x  6  fl,  t  >  0 

,  i 
1=  1 

where  A  =  X{t)  is  a  Lagrange  multiplier  associated  to  the  constraint  (2.5c),  A 
is  a  bounded  linear  operator  from  L^(Q)  into  and  A*  denotes  its  adjoint. 

Finally,  a,(p)  =  ^o(p)  (i  =  1,2)  where  a  is  smooth  (say,  with  bounded 

derivatives  up  to  order  2),  a  is  spherically  symmetric  (rotational  invariance).  Of 
course,  the  model  introduced  in  this  work  corresponds  to  a(p)  =  |p|.  However, 
this  choice  induces  such  singularities  that  a  mathematical  analysis  does  not  seem 
to  be  possible.  This  is  why  we  shall  study  some  model  cases  involving  slightly 
regularized  \’ariants  of  this  choice.  In  fact,  numerically,  scaling  (2.9a)  also  induces 
some  regularizations  of  a  quite  similar  to  the  ones  we  make  below  -  and,  thus,  our 
analysis  covers  situations  that  are  quite  realistic  from  the  numerical  viewpoint. 

We  shall  therefore  assume  that  there  exists  £  (0, 1)  such  that 
(5.3)  <(iiii,)foraJlp€R2 

in  the  sense  of  symmetric  matrices. 

If  we  go  back  to  our  real  choice  of  a,  namely  a(p)  =  |p|,  we  see  that  (5.3)  does 
not  hold  for  p  near  zero  and  for  |p|  — >  oo.  The  singularity  at  p  =  0  induces 
mathematical  and  numerical  difficulties.  In  practice  we  truncate  near  p  =  0. 
The  assumption  for  p  large  can  be  relaxed  by  proving  some  upper  bounds  with  a 
rather  technical  argument  (contained  in  the  proof  below).  We  prefer  to  skip  this 
technical  argument  in  order  to  avoid  confusing  the  main  issue  which  concern  the 
imposition  of  contraints. 

Finally,  we  observe  that  enforcing  (2.5c)  while  scaling  (5.2a)  amounts  to  requiring 
that 

(5.2b)  A  =  >l[u]  A*{Au  -  uq  j<ij)  \\A'‘{Au  -  \\Au  -  uo||i*(n)  = 

where  ^[u]  =  ^(a,(Vu)). 

We  also  prescribe  an  initial  condHion 

(5.2c)  u]t=ro  =  in  fl 
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zind  boundary  condition 


(.'i.2d) 


du 

dn 


=  0  on  dQ. 


Of  course,  we  need  to  assume  that  u°  satisfies 


(5.4) 


IIAu"  -  Uo||LJ(n)  =  <T- 


Remark  5.2;  The  existence  of  some  u°  satisfying  (5.4)  is  not  obvious  zmd  de¬ 
pends  very  much  on  the  properties  of  A.  This  existence  is  certainly  ensured  by  the 
following  assumption 

(o.5)  i?(>l)(range  of  A)  is  dense  in 

Indeed,  this  condition  implies  that  there  exist  G  (or  as  smooth  as  we 

wish  by  density  and  continuity  of  >1)  satisfying 

(5.6)  ||Au^  -  uoWlhh)  <  \\Au‘^  -  uo|U2(ft)  >  a. 

And  it  is  enough  to  tadce  -f  (1  —  6)u^  for  some  convenient  0  €  (0, 1).  □ 

In  order  to  illustrate  clearly  the  mathematical  difficulties  and  results  associated 
with  the  system  (5.2a-d)  we  begin  with  the  model  case  where  a{p)  =  so  that 
ai(Vu)  =  Vu  and  A(u]  reduces  to  the  classical  lineeir  operator  A[u]  =  Au  (the 
Laplace  operator).  In  that  case,  we  can  prove  the 

Theorem  5.1:  Let  A  satisfy  (5.5),  let  u°  €  /f^(S7)  satisfy  (5.4),  let  uo  €  L^(fl)  and 
let  a(p)  =  |1pP  on  R^.  Then,  there  exists  a  unique  solution  u  o/(5.2a-d)  satisfying: 
u  €  Ci[0,ooy,H^n))nL^i0,T-H^{Q)),  Ut  €  l2(0,T;L2(n)),  A  G  l'"(0,r)  for  all 
T  G  (0,  oo). 

Remark  5.3:  If  A  is  a  convolution  operator  i.e.  Au  =  k  *u  then  (5.5)  holds  if  and 
only  if  the  Fourier  transform  k  of  k  satisfies 

(5.7)  meas  {(  G  RVMO  =  0}  =  0. 


Remark  5.4:  If  A  and  A*  map  boundedly  into  for  2  <  p  <  oo  then 

the  proof  below  also  show  that  u  G  C([0,  oo);  IF^’'’(fI))  D  L^{0,T’,  ut  G 

X»’(0,T;XJ’(fI)),  A  G  1^(0, T)  for  all  T  G  (0,oo)  if  u°  G  ^^'''’(n),  uo  € 
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Remark  5.5:  The  proof  also  applies  to  different  type  of  constraints  that  can  even 
be  nonquadratic  constraints.  Let  us  only  mention  a  few  possibilities  for  which  the 
same  result  as  the  one  above  holds.  In  the  case  of  multiplicative  noise  as  in  the 
previous  section,  we  can  replace  (5.2b)  by 

(5.8)  J  dx  =  (7^  >  0 

assuming  for  instance  that  6  £°°(fl).  Also,  we  might  want  to  enforce  local 
constraints  on  a  finite  peu'tition  (or  subpartition)  of  fl,  that  in  practice  can  be 
obtained  by  a  segmentation  algorithm.  In  that  case,  we  consider  wi , . . .  ,  u}m(^  ^  1) 
measurable  sets  in  such  that  measfwi  0  u>j]  =  0  for  all  1  <  i  ^  j  <  m  and  we 
replace  (5.2b) 

(5.9)  /  |Au  —  =  <7^  >  0. 


Then,  Theorem  (5.1)  still  holds  for  the  corresponding  condition  equation  that  in¬ 
volves  now  m  different  Lagrange  multipliers  A,  =  A[u]  A*{Au  —  uo)dx.  □ 

Proof  of  Theorem  5.1: 

Step  1:  General  a  priori  estimates. 

Here,  we  list  some  general  consquences  of  the  fact  that  the  evolution  equation  we 
are  considering  is  a  gradient  flow  (of  a  constreiined  functional).  Indeed,  multiplying 
(5.2a)  by  ut  and  using  (5.2b),  we  deduce 

(5.10)  f  \ut\^dx  +  -Y  f  a(Vu)dx  =  0  for  t  >  0. 

Ju  dt  Jq 

Hence,  Ut  is  bounded  in  L^(0,  oo;  Z,^(I7))  and  Vu  is  bounded  in  L°°(0,  cc;  Z-^(f2)). 
In  pEirticular,  u  =  fluids  -|-  u°  is  bounded  in  for  all  T  €  (0, oo) 

and  thus  u  is  bounded  in  C([0,  T]; /f*(n))  for  all  T  €  (0,  oo). 

Step  2:  A  lower  bound 

We  want  to  show  that  \\A*{Au  —  uo)||z,2(n)  is  bounded  from  below  uniformly  on 
(0,T]  (for  all  T  e  (0,  oo))  and  that  the  lower  bound  depends  only  on  T  and  on  the 
norm  of  u“.  Indeed,  if  this  were  not  the  case,  in  view  of  the  estimates  shown  in 
Step  1  and  in  view  of  (5.2b),  this  would  yield  the  existence  of  a  sequence  {uy}j>i 
such  that  Uj  is  bounded  in  and 


(5.11)  II  Auj  -  uo||LJ(n)  =  \\A*{Auj  -  uo)||L2(n)  y  0. 
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Without  loss  of  generality,  we  may  assume,  extracting  a  subsequence  if  necessary, 
that  Uj  converges  weakly  in  H^(Q)  to  some  u  and  thus  by  the  Rellich-Kondrakov 
theorem,  Uj  converges  strongly  in  to  some  u.  Since  A  and  A*  are  bounded 

from  into  (5.11)  then  implies 


(5.12)  l|Au  —  uolli,*(n)  =  A*{Au  —  uo)  =  0. 

In  other  words,  Au  —  uq  belongs  to  the  kernel  of  A*.  But  (5.5)  implies  that  this 
kernel  is  trivial  (reduces  to  {0})  therefore  Au  =  uq  and  we  reach  a  contradiction 
with  the  first  statement  in  (5.12). 

Step  S:  estimates. 

We  multiply  (5.2a)  by  —Au  and  we  find 

(5.13) 

~  J^^\Vu\^dxA  J^{-Aufdx  =  {J^{-Au)[A%Au-uo)]dxf{\\A*{Au-uo)UHU))-^. 

We  then  fix  T  €  (0,  oo),  use  elliptic  regularity  and  Steps  1  and  2  to  deduce 

T 

(5.14)  l|u|li=(o,T,H>(!l))  £  Co(l  +  /  1  /  (-A«)(.4'(.4u  -  u,)]dxl'‘dt) 

Jq  Jn 

where  Co  depends  only  on  T  and  on  bounds  of 

Next,  we  observe  that  since  {u{t)/t  6  [0,  Tj}  is  bounded  in  ff^(fi),  by  Rellich- 
Kondrakov  theorem,  {u{t)/t  e  [0, T]}  is  relatively  compact  in  £^(0)  and  thus  since 
A  and  A*  are  bounded  from  L^{Q,)  into  I<^(ft),  {A*(Au(<)  —  uo)/t  6  [0,T]}  is 
relatively  compact  in  L^(f2).  This  implies  that  we  can  decompose  A*{Au{t)  —  uq) 
as  follows  for  all  e  >  0 

(5.15a)  A*{Auit)-uo)  =  f{t)  +  git)  V<g[0,T] 

(5.15b)  \\m\\LHn)<s,  ll^(0ll//-(n)  <  Vfe[0,r] 


for  some  C(e)  that  depends  only  on  e,T  emd  bounds  of  u“. 
Therefore,  we  have  for  all  t  G  [0,r] 


(-Au){A*(Au  -  uo)}dx|  <  e|lul|//2(n)  +  1  / 

Jn 

S  £||«ll//»(n)  +  I  f 

Jn 


{-Au)g{t)dx\ 
Vu  ■  V9(t)d'x|. 
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Since  u  satisfies  (5. 2d),  and  thus  finally 

I  /  (-Au){>l*(^u  -  uo)}dil  <  e||ull//*(n)  +  C'(£)||u||j:/»(n)- 

Ju 

Hence,  if  we  input  this  bound  in  (5.14)  and  use  the  bound  shown  in  Step  1,  we 
deduce 

ll“lliHo,T;H*(n))  ^  +  2C'oellu||i2(o,r;H*(n)) 

and  we  conclude  by  choosing  £  =  . 

Step  4-  Uniqueness. 

We  consider  two  solutions  u,v  of  (5.2a-d)  and  denote  by  A,/i  the  corresponding 
Lagrange  multipliers.  Obviously,  we  have,  for  u  —  v  =  w 

Wt  —  Aw  =  A  A*Aw  =  (A  —  pl)A*{Av  —  uo). 

Multiplying  this  equation  by  w  and  — Aa>,  integrating  by  parts  and  summing  up, 
we  find  easily  for  all  To(0,  oo) 

(5.16) 

il“’(‘)llHHn)  +  ll“’lli>(0,iiH>(n))  f  X{s)ds  j  {A‘Xu))(-Ati>  +  u')iii|+ 

Jo  Jn 

+  f  jA  —  ^|dsl  j  [A*(Au  —  uo)}(— Alt)  +  iy)<iil}  for  alH  e  [0,T) 

Jo  Ja 

for  some  positive  constant  Ci  depending  only  on  T. 

Using  the  same  argument  as  in  Step  3,  we  deduce  the  following  bounds  for  eJI 
£  >  0 

(5.17a)  I  / (yl*Au))(-Au)  +  u))di|  <£|lu)l|H2(n)||i"IU>(n)  +  C'(£)||u)||^i(n) 

Jci 

(5.17b)  )  [  {A*{Av  -  uo)}(-Au)  +  w)dx\  <  e)!u;))//2(n)  +  C(£)))u))|//i(n) 

Jsi 

where  C(£)  denotes  various  positive  constants  depending  only  on  £  and  T.  Inserting 
these  bounds  in  (5.16)  we  find 

Jo 

+  Cl  j  lA  - /xl{£l|it)ll//2(n)  +  C(e)llu;||//i(n)}ds. 

Jo 
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Using  the  Cauchy-Schwarz  inequzility,  we  deduce  easily 
(5.18) 

ll^(OllHi(n)  +  ll“^lli2(o,<;i/*(n))  -  ^ll^llL»(o,f;//*(n))  +C'(£  /  a(^)||it^('S)llHi(n)ds+ 

Jo 


+Ci  /  lA  - /il{£|lu;||//j(n)  +  C(£)||u)ll//i(n)}dfi 
Jo 


where  a  >  0,  a  €  £^(0,  T)  (for  all  T  €  (0,  oo)). 

(5-19) 

Next,  we  estimate  \  —  fi.  In  view  of  Step  2,  we  have 

|A  -  |xl  <  C2[l|ii)||L2(n)|  /  (-Au){i4*Au  -  uo)}dx|  +  I  /  Awdx\+ 

Jn  Jq 

+  1  /  (— Aiy){A*>lu  —  uo)}da:| 

Ja 

for  some  C2  >  0  which  depends  only  on  T.  But  this  yields  immediately 
|A  -  )u|  <  61|u;|(£,s(n)  4-  C3||tx;|(//2(n),  for  some  b  €  L^(0,T). 

Going  back  to  (5.18),  we  obtmn  for  all  <  e  [0,T]. 

+  Il^lli2(0, <;//*(«))  ^  ^(1  +  ^1^3)1|u)|||2(0  + 

C(£)  /  o(s)l|u;(s)||]^.(njds +  CiC'3C'(e)  /  ||u;ll//i(n)||u;|j//2(n)ds+ 

Jo  Jo 

+  Ci£  f  6('S)ll^llL2(n)ll«'llH2(n)<i'S  +  CiC(£)  f  fe(^)ll“’llL2(n)ll’^llH>(n)‘^’®- 
Jo  Jo 

Using  the  Cauchy- Schwaxz  inequality,  this  yields: 

ll^(0llH>(n)  +  li*"llL2(o,«;H2(n))  <€^(2  +  ^i<^3)||w||i2(o,t;//2(n)+ 

/  c('S)|(u)(s)||^,(jj)‘^^ 

Jo 

where  c  >  0,  c  £  £*(0,  T). 
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We  then  choose  £  =  (2  +  C1C3)  ^  and  we  conclude  using  Gronwall’s  inequality. 

Conclusion:  We  conclude  the  proof  of  Theorem  5.1  here  since  only  the  existence 
part  has  not  been  completed.  But  this  part  is  a  straightforward  consequence  of 
solving  approximate  problems,  proving  the  same  a  priori  bounds  uniformly  for  the 
approximate  solution  and  passing  to  the  limit  using  these  bounds.  We  do  not  want 
to  give  all  the  details  of  such  a  tedious  argument.  Let  us  only  mention  a  few  possible 
approximations  like  a  penalty  method  (penalizing  the  constraint),  implicit  time 
discretization  (solving  each  stationary  problem,  for  each  iteration,  by  a  minimization 
problem  similar  to  the  ones  solved  in  Proposition  5.1),  or  splitting  methods  similar 
to  the  numerical  method  presented  in  the  following  section  (where  we  solve  first  the 
equation  without  constraints  on  a  time  interval  of  length  A<  and  we  then  project 
back  to  the  obtained  solution  on  the  constraints  manifold  by  a  simple  affine  rule). 

For  all  these  approximation  methods,  one  can  adapt  the  a  priori  estimates  shown 
above.  But  we  certainly  do  not  want  to  do  so  here  in  order  to  avoid  confusing  the 
main  issues  in  this  paper.  □ 

We  now  turn  to  a  nonlinear  equation  (5.2a)  where  a(p)  satisfies  the  conditions 
mentioned  in  the  beginning  of  this  section  and  (5.3)  in  particular. 

Theorem  5.2.  Let  a{p)  satisfy  (5.3),  let  u®  €  satisfy  (5.4),  let  uq  £  L^(fl) 

and  let  A  satisfy  (5.5).  Then, 

i)  there  exists  a  solution  u  o/(5.2a-d)  satisfying:  u  €  C([0, 00);  i?^(fl))njL^(0,  T;  if^(n)), 
ut  €  1^(0,  T;  L2(f2)),  A  €  T)  for  all  T  €  (0, 00). 

ii)  If  Uo  £  H^{Q,)  and  A,  A*  are  bounded  from  H^{Q)  into  then  the 

solution  is  unique  and  A  £  L°°{0,T)  for  all  T  £  (0, 00). 

Remark  5.6.  The  analogues  of  Remarks  5.4  and  (5.5)  hold  here.  In  particular, 
using  this  extra  regularity,  one  can  show  the  uniqueness  of  solutions  by  an  argument 
quite  similar  to  the  one  given  in  the  proof  of  Theorem  3.1  euid  which  does  not  use 
a  regularity  assumption  on  uq  like  in  part  ii)  of  Theorem  3.2  above.  However,  this 
argument  relies  upon  a  regularity  result  which  is  too  technical  to  be  detailed  here. 

Remark  5.7.  If  fl  =  and  A  is  a  convolution  operator  then  A,  A*  are  bounded 
operators  from  into  since  they  are  bounded  from  L^{Q)  into  L^(fl) 

and  they  commute  with  differentiation. 

Proof  of  Theorem  5.2.  We  only  explain  the  modifications  that  have  to  be  made 
in  the  proof  of  Theorem  5.1.  In  particular.  Steps  1  and  2  are  identical.  However, 

Step  3  has  to  be  modified  substantially.  The  final  result  being  the  same,  these  facts 
and  the  uniqueness  argument  shown  below  allow  us  to  complete  the  proof  of  part 

i)- 

The  L^{Hl)  bound  follows  from  multiplying  (5.2a)  by  —  Au  and  making  some 
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observations.  First  of  all,  >l[u]  =  dx^dxj  ''^^ere  Oij  =  and 

thus  l>t[u]|  <  C\D^u\.  Next,  we  recall  an  adaptation  of  a  famous  inequality  due  to 
H.O.  Cordes  [16]  (see  also  A. I.  Koselev  [17])  shown  in  P.L.  Lions  [18]  in  the  case 
of  Neumann  boundary  conditions.  There  exist  a  >  0,  C  >  0  such  that  for  edl 
u  €  satisfying  (3. Id) 

(5.19)  /  ^[u]  Audi  >  Q;(|u||^,(n)  -  C’||ul|^2(n). 

Jn 

This  inequality  allows  us  to  adapt  easily  the  rest  of  the  proof  made  in  Step  3. 

We  conclude  with  the  proof  of  the  uniqueness  statement  (part  ii)  above.  Let  u,  v 
be  two  solutions  of  (5.2a-d)  and  let  A,/i  be  the  corresponding  Lagrange  multipliers. 
We  denote  hy  w  =  u  —  v  and  multiply  by  w  the  equation  satisfied  by  w.  We  then 
find  in  view  of  (5.3) 

(5.20)  f  w^dx  +  u  f  |Vt£jpdi  <  C{|A|  f  w^dx  +  [A  —  ^|(  /  w^dx)^}. 

2  dt  Jq  Jq  Jq  Jq 

Next,  we  observe  that  A  can  be  written  as 

(5.21)  A  =  -{/  ^ai(Vu)— {A*(Au-uo)}dx}llA*(Au-uo)|li?(n) 

JCl 

and  a  simiW  expression  holds  for  n  with  u  replaced  by  v. 

From  these  expressions  we  deduce  using  the  assumptions  made  about  a.  A,  A* 
and  uo 

(5.22)  |A-Ml<qHI//>(n) 

(recall  that  u,v  are  bounded  in  L°°{Q,T^H^{U))  for  all  T  €  (0,  oo)). 

Inserting  this  bound  in  (5.20)  we  finally  deduce 

2^ll“^lli*(n)  -  ^{|A|||«’lli2(n)  +  li“^llH>(n)!l«^llL2(n)- 

Hence,  we  have  for  all  t  G  [0,  T],  by  the  Cauchy-Schwarz  inequality 
(5-23)  +  lADllii’lll^cn) 
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and  the  uniqueness  follows  from  Gronwall’s  inequality. 

Let  us  finally  observe  that  the  fact  that  A  6  L°°(0,  T)  (for  all  T  E  (0,  oo)) 
is  straightforward  in  view  of  (5.21)  since  u  is  bounded  in  and  \\A*{Au  — 

uo)l|L2(n'>  is  bounded  from  below  (Step  2).  □ 

6.  Results. 

Page  PI  shows  an  extremely  noisy  image  (Picture  2)  with  various  standard  meth¬ 
ods  which  we  compare  with  our  TV  based  method.  Page  P2  shows  the  result  of 
edge  detection  applied  to  the  images  on  Pi.  Notice  picture  P12  -  the  TV  denoised 
image  admits  excellent  edge  detection. 

Pages  P3  and  P4  show  restoration  when  the  noise  is  multiplicative  as  described  in 
section  4.  The  only  information  we  use  concerns  the  mean  and  standard  deviation 
of  the  noise.  No  other  method  can  handle  these  situations. 

Page  P5  shows  the  “tank”  restoration  and  pages  P6-7  show  the  “chemical  plant” 
restoration. 

On  page  P8,  we  begin  our  deblurring/ denoising  demonstration.  We  use  various 
types  of  blur  with  various  amounts  of  additive  white  noise.  Our  restorations  are 
shown  on  the  right.  On  page  P9  we  show  the  results  of  other  methods  applied  to 
the  motion  blurred  image,  on  page  PIO  we  do  it  for  the  Gaussian  blur,  euid  on  page 
PI  1  we  do  it  for  the  out  of  focus  blur. 

We  believe  that  our  restoration  technique,  supported  by  this  contract,  represents 
a  break  through  in  the  area. 
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